This analysis
assesses how felzartamab therapy affects molecular sores measured with
gene expression data from biopsies from kidney allografts taken from
patients with antibody mediated rejection. The assumptions of normality
and heterogeneity of variance differ across molecular scores, thus we
avail to alinged-rank transformation to carry out a non-parametric
ANOVA. This is powerful in this context because it allows us to test for
the interactive effect of treatment (i.e., placebo vs felzartamab) and
time (i.e., change in scores from biopsy to biopsy). Aligned-rank
transformation also allows for us to specify a mixed-model to handle
repeated measurements (i.e., biopsies) across patients.
References:
• Wobbrock J, F.L., Gergle D, Higgins J The Aligned Rank Transform for Nonparametric Factorial Analyses Using Only ANOVA Procedures. In Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ’11) 143–146. doi:10.1145/1978942.1978963, https://depts.washington.edu/acelab/proj/art/.(2011).
• Kay M, E.L., Higgins J, Wobbrock J. ARTool: Aligned Rank Transform for Nonparametric Factorial ANOVAs. doi:10.5281/zenodo.594511, R package version 0.11.1, https://github.com/mjskay/ARTool. (2021).
• Elkin, L.A., Kay, M., Higgins, J.J. & Wobbrock, J.O. An Aligned Rank Transform Procedure for Multifactor Contrast Tests. in The 34th Annual ACM Symposium on User Interface Software and Technology 754–768 (Association for Computing Machinery, Virtual Event, USA, 2021).
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# Load necessary libraries
library(tidyverse)
library(broom)
library(rstatix)
library(flextable)
library(officer)
library(emmeans)
library(multcomp)
library(ARTool)
library(rcompanion)
library(Biobase)
library(ggpubr)
library(patchwork)
library(ggprism)
# Custom operators
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# Suppress dplyr reframe info
options(dplyr.reframe.inform = FALSE)
# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/results/flextable_artANOVA.RData")
# load plot data
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_artANOVA.RData")
# Seed for reproducibility
set.seed(42)
Now we
define the categories of molecular scores for easier organization
downstream.
# Feature categories
vars_cfDNA <- c("cfDNA_cpml")
vars_abmr <- c("ABMRpm", "ggt0", "ptcgt0", "DSAST", "NKB")
vars_tcmr <- c("TCMRt", "tgt1", "igt1", "QCAT", "TCB")
vars_injury <- c("IRRAT30", "IRITD3", "IRITD5", "cigt1", "ctgt1")
vars_parenchyma <- c("KT1", "KT2")
vars_macrophage <- c("AMAT1", "QCMAT")
# DEFINE VARIABLES TO ASSESS ####
vars <- c(vars_cfDNA, vars_abmr, vars_tcmr, vars_macrophage, vars_injury, vars_parenchyma)
Now we
process the data so we can iteratively carry out the analyses on each
score.
# DEFINE THE DATA ####
data <- set %>%
pData() %>%
dplyr::select(Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, all_of(vars)) %>%
dplyr::filter(Patient %nin% c(15, 18)) %>%
left_join(., summarise(., sample_pairs = n(), .by = Felzartamab_Followup), by = "Felzartamab_Followup") %>%
relocate(sample_pairs, .after = Felzartamab_Followup) %>%
arrange(Felzartamab, Patient, Group)
# WRANGLE THE PHENOTYPE DATA ####
data_00 <- data %>%
expand_grid(category = c("cfDNA", "ABMR", "TCMR", "macrophage", "injury", "parenchyma")) %>%
nest(.by = category) %>%
mutate(
features = map(
category,
function(category) {
if (category == "cfDNA") {
vars_cfDNA
} else if (category == "ABMR") {
vars_abmr
} else if (category == "TCMR") {
vars_tcmr
} else if (category == "macrophage") {
vars_macrophage
} else if (category == "injury") {
vars_injury
} else if (category == "parenchyma") {
vars_parenchyma
}
}
),
data = pmap(
list(features, data),
function(features, data) {
data %>%
dplyr::select(
Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, sample_pairs,
all_of(features)
)
}
)
)
# WRANGLE THE DATA FOR UNIVARIATE TESTS ####
data_01 <- data_00 %>%
mutate(
data_univariate = pmap(
list(data, features),
function(data, features) {
data %>%
pivot_longer(
cols = all_of(features),
names_to = "variable",
values_to = "value"
) %>%
nest(.by = variable) %>%
mutate(
variable = variable %>%
factor(
levels = vars
),
annotation = case_when(
variable %in% vars_cfDNA ~ "cfDNA",
variable %in% vars_tcmr ~ "TCMR-related",
variable %in% vars_abmr ~ "ABMR-related",
variable %in% vars_macrophage ~ "macrophage-related",
variable %in% vars_injury ~ "injury-related",
variable %in% vars_parenchyma ~ "parenchyma-related",
TRUE ~ " "
) %>%
factor(
levels = c(
"cfDNA",
"ABMR-related",
"TCMR-related",
"macrophage-related",
"injury-related",
"parenchyma-related",
"archetypes",
"rejection PC",
"injury PC"
)
),
score = case_when(
variable == "cfDNA_cpml" ~ "Donor-derived cell-free DNA (dd-cfDNA, cp/mL)",
variable == "TCMRt" ~ "TCMR classifier (TCMRProb)",
variable == "TCB" ~ "T cell burden (TCB)",
variable == "tgt1" ~ "Tubulitis classifier (t>1Prob)",
variable == "igt1" ~ "Interstitial infiltrate classifier (i>1Prob)",
variable == "QCAT" ~ "Cytotoxic T cell-associated (QCAT)",
variable == "ABMRpm" ~ "ABMR classifier (ABMRProb)",
variable == "DSAST" ~ "DSA-selective transcripts (DSAST)",
variable == "NKB" ~ "NK cell burden (NKB)",
variable == "ggt0" ~ "Glomerulitis classifier (g>0Prob)",
variable == "ptcgt0" ~ "Peritubular capillaritis classifier (ptc>0Prob)",
variable == "AMAT1" ~ "Alternatively activated macrophage (AMAT1)",
variable == "QCMAT" ~ "Constitutive macrophage (QCMAT)",
variable == "IRITD3" ~ "Injury-repair induced, day 3 (IRITD3)",
variable == "IRITD5" ~ "Injury-repair induced, day 5 (IRITD5)",
variable == "IRRAT30" ~ "Injury-repair associated (IRRAT30)",
variable == "cigt1" ~ "Fibrosis classifier (ci>1Prob)",
variable == "ctgt1" ~ "Atrophy classifier (ct>1Prob)",
variable == "KT1" ~ "Kidney parenchymal (KT1)",
variable == "KT2" ~ "Kidney parenchymal - no solute carriers (KT2)"
), .before = 1
) %>%
arrange(annotation, variable)
}
)
) %>%
dplyr::select(category, data_univariate) %>%
unnest(data_univariate) %>%
mutate(n_cat = score %>% unique() %>% length(), .by = category, .after = variable)
This is
what the data look like after we have organized them. Each row of the
dataframe contains the data for each score, the category it belongs to,
its category annotation, its long-format name, and the name of the
variable in the data.
data_01 %>% print(n="all")
## # A tibble: 20 x 6
## category annotation score variable n_cat data
## <chr> <fct> <chr> <fct> <int> <list>
## 1 cfDNA cfDNA Donor-derive~ cfDNA_c~ 1 <tibble>
## 2 ABMR ABMR-related ABMR classif~ ABMRpm 5 <tibble>
## 3 ABMR ABMR-related Glomerulitis~ ggt0 5 <tibble>
## 4 ABMR ABMR-related Peritubular ~ ptcgt0 5 <tibble>
## 5 ABMR ABMR-related DSA-selectiv~ DSAST 5 <tibble>
## 6 ABMR ABMR-related NK cell burd~ NKB 5 <tibble>
## 7 TCMR TCMR-related TCMR classif~ TCMRt 5 <tibble>
## 8 TCMR TCMR-related Tubulitis cl~ tgt1 5 <tibble>
## 9 TCMR TCMR-related Interstitial~ igt1 5 <tibble>
## 10 TCMR TCMR-related Cytotoxic T ~ QCAT 5 <tibble>
## 11 TCMR TCMR-related T cell burde~ TCB 5 <tibble>
## 12 macrophage macrophage-related Alternativel~ AMAT1 2 <tibble>
## 13 macrophage macrophage-related Constitutive~ QCMAT 2 <tibble>
## 14 injury injury-related Injury-repai~ IRRAT30 5 <tibble>
## 15 injury injury-related Injury-repai~ IRITD3 5 <tibble>
## 16 injury injury-related Injury-repai~ IRITD5 5 <tibble>
## 17 injury injury-related Fibrosis cla~ cigt1 5 <tibble>
## 18 injury injury-related Atrophy clas~ ctgt1 5 <tibble>
## 19 parenchyma parenchyma-related Kidney paren~ KT1 2 <tibble>
## 20 parenchyma parenchyma-related Kidney paren~ KT2 2 <tibble>
At this
stage we can summarize the median values for each score by treatment and
follow-up.
(Note that the median scores are used for general interpretation of
treatment effects because aligned-ranks, the data being assessed for
inference regarding the effects of treatment, are not easily
interpretable.)
# UNIVARIATE MEDIANS ####
data_02 <- data_01 %>%
mutate(
medians = map(
data,
function(data) {
data %>%
reframe(
median = value %>% median(),
IQR = value %>% IQR(),
sample_pairs = n(),
.by = c(Followup, Felzartamab, Felzartamab_Followup)
) %>%
arrange(Felzartamab_Followup)
}
),
medians_delta = map(
medians,
function(medians) {
medians %>%
reframe(
median_delta = combn(median, 2, diff) %>% as.numeric(),
IQR_delta = combn(IQR, 2, mean) %>% as.numeric(),
sample_pairs = sample_pairs,
.by = c(Felzartamab)
) %>%
mutate(
Followup_pairwise = rep(c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52"), 2) %>%
factor(levels = c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52")),
.before = sample_pairs
) %>%
mutate(
median_delta_delta = combn(median_delta, 2, diff) %>% as.numeric(),
IQR_delta_delta = combn(IQR_delta, 2, mean) %>% as.numeric(),
.by = c(Followup_pairwise),
.before = sample_pairs
)
}
)
)
Let’s
take a look at what the median summaries look like for the ABMRprob
score.
data_02 %>%
dplyr::filter(variable == "ABMRpm") %>%
pull(medians) %>%
pluck(1) %>%
flextable::flextable()
Followup | Felzartamab | Felzartamab_Followup | median | IQR | sample_pairs |
|---|---|---|---|---|---|
Baseline | Placebo | Baseline_Placebo | 0.6695 | 0.35525 | 10 |
Week24 | Placebo | Week24_Placebo | 0.7740 | 0.39025 | 10 |
Week52 | Placebo | Week52_Placebo | 0.5185 | 0.49825 | 10 |
Baseline | Felzartamab | Baseline_Felzartamab | 0.7360 | 0.30775 | 10 |
Week24 | Felzartamab | Week24_Felzartamab | 0.1675 | 0.35700 | 10 |
Week52 | Felzartamab | Week52_Felzartamab | 0.7025 | 0.53875 | 10 |
Now we
can look at the change in median ABMRprob scores in the placebo and
felzartamab arms, as well as the differnce in those medians (i.e., ΔΔ
median - this is our summarized treatment effect).
data_02 %>%
dplyr::filter(variable == "ABMRpm") %>%
pull(medians_delta) %>%
pluck(1) %>%
flextable::flextable()
Felzartamab | median_delta | IQR_delta | Followup_pairwise | median_delta_delta | IQR_delta_delta | sample_pairs |
|---|---|---|---|---|---|---|
Placebo | 0.1045 | 0.372750 | Baseline - Week24 | -0.6730 | 0.3525625 | 10 |
Placebo | -0.1510 | 0.426750 | Baseline - Week52 | 0.1175 | 0.4250000 | 10 |
Placebo | -0.2555 | 0.444250 | Week24 - Week52 | 0.7905 | 0.4460625 | 10 |
Felzartamab | -0.5685 | 0.332375 | Baseline - Week24 | -0.6730 | 0.3525625 | 10 |
Felzartamab | -0.0335 | 0.423250 | Baseline - Week52 | 0.1175 | 0.4250000 | 10 |
Felzartamab | 0.5350 | 0.447875 | Week24 - Week52 | 0.7905 | 0.4460625 | 10 |
Now we
carry out the analysis. This includes testing the specific contrasts for
the interaction of treatment and time (i.e., our effect of
treatment).
# UNIVARIATE NONPARAMETRIC TESTS ####
artANOVA <- data_02 %>%
mutate(
art = map(data, art, formula = value ~ Followup * Felzartamab + (1 | Patient)),
art_aov = map(art, anova, type = "II"),
art_aov_tidy = map(art_aov, tidy) %>% suppressWarnings(),
art_aov_contrast = map(
art,
art.con,
formula = "Followup:Felzartamab", adjust = "fdr", method = "pairwise", interaction = TRUE, response = "art"
),
art_aov_contrast_tidy = map(art_aov_contrast, tidy),
art_aov_contrast_cld = map(
art_aov_contrast_tidy,
function(art_aov_contrast_tidy) {
art_aov_contrast_tidy %>%
as.data.frame() %>%
cldList(adj.p.value ~ Followup:Felzartamab, data = .)
}
)
)
# CREATE FLEXTABLE SUMMARY OF ART MODELS ####
title_art <- paste("Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients")
artANOVA %>%
dplyr::select(annotation, n_cat, score, art_aov) %>%
unnest(everything()) %>%
dplyr::rename(p.value = `Pr(>F)`) %>%
mutate(FDR = p.value %>% p.adjust(method = "fdr"), .by = annotation) %>%
dplyr::select(annotation, score, Term, `F`, p.value, FDR) %>%
flextable::flextable() %>%
flextable::add_header_row(values = rep(title_art, ncol_keys(.))) %>%
flextable::merge_h(part = "header") %>%
flextable::merge_v(j = 1:2) %>%
flextable::fontsize(size = 8, part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::bg(i = ~ FDR < 0.05 & Term == "Followup:Felzartamab", j = 3:6, bg = "#fbff00") %>%
flextable::colformat_double(j = 2:3, digits = 2) %>%
flextable::colformat_double(j = 4:ncol_keys(.), digits = 3) %>%
flextable::border_remove() %>%
flextable::bold(part = "header") %>%
flextable::padding(padding = 0, part = "all") %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::autofit()
Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients | |||||
|---|---|---|---|---|---|
annotation | score | Term | F | p.value | FDR |
cfDNA | Donor-derived cell-free DNA (dd-cfDNA, cp/mL) | Followup | 4.418 | 0.019 | 0.029 |
Felzartamab | 1.598 | 0.222 | 0.222 | ||
Followup:Felzartamab | 10.396 | 0.000 | 0.001 | ||
ABMR-related | ABMR classifier (ABMRProb) | Followup | 6.659 | 0.003 | 0.007 |
Felzartamab | 1.620 | 0.219 | 0.235 | ||
Followup:Felzartamab | 9.605 | 0.000 | 0.001 | ||
Glomerulitis classifier (g>0Prob) | Followup | 11.768 | 0.000 | 0.001 | |
Felzartamab | 1.511 | 0.235 | 0.235 | ||
Followup:Felzartamab | 10.743 | 0.000 | 0.001 | ||
Peritubular capillaritis classifier (ptc>0Prob) | Followup | 14.188 | 0.000 | 0.000 | |
Felzartamab | 3.178 | 0.092 | 0.125 | ||
Followup:Felzartamab | 7.749 | 0.002 | 0.004 | ||
DSA-selective transcripts (DSAST) | Followup | 10.508 | 0.000 | 0.001 | |
Felzartamab | 1.747 | 0.203 | 0.234 | ||
Followup:Felzartamab | 4.179 | 0.023 | 0.035 | ||
NK cell burden (NKB) | Followup | 6.512 | 0.004 | 0.007 | |
Felzartamab | 1.956 | 0.179 | 0.224 | ||
Followup:Felzartamab | 5.114 | 0.011 | 0.018 | ||
TCMR-related | TCMR classifier (TCMRProb) | Followup | 3.555 | 0.039 | 0.165 |
Felzartamab | 1.417 | 0.249 | 0.468 | ||
Followup:Felzartamab | 1.124 | 0.336 | 0.504 | ||
Tubulitis classifier (t>1Prob) | Followup | 1.016 | 0.372 | 0.506 | |
Felzartamab | 1.853 | 0.190 | 0.408 | ||
Followup:Felzartamab | 0.065 | 0.937 | 0.968 | ||
Interstitial infiltrate classifier (i>1Prob) | Followup | 3.414 | 0.044 | 0.165 | |
Felzartamab | 0.997 | 0.331 | 0.504 | ||
Followup:Felzartamab | 0.032 | 0.968 | 0.968 | ||
Cytotoxic T cell-associated (QCAT) | Followup | 5.693 | 0.007 | 0.107 | |
Felzartamab | 4.118 | 0.057 | 0.172 | ||
Followup:Felzartamab | 0.927 | 0.405 | 0.506 | ||
T cell burden (TCB) | Followup | 4.360 | 0.020 | 0.151 | |
Felzartamab | 2.868 | 0.108 | 0.269 | ||
Followup:Felzartamab | 0.301 | 0.742 | 0.856 | ||
macrophage-related | Alternatively activated macrophage (AMAT1) | Followup | 1.549 | 0.226 | 0.304 |
Felzartamab | 1.047 | 0.320 | 0.320 | ||
Followup:Felzartamab | 2.406 | 0.105 | 0.209 | ||
Constitutive macrophage (QCMAT) | Followup | 3.378 | 0.045 | 0.209 | |
Felzartamab | 3.100 | 0.095 | 0.209 | ||
Followup:Felzartamab | 1.426 | 0.253 | 0.304 | ||
injury-related | Injury-repair associated (IRRAT30) | Followup | 0.055 | 0.946 | 0.946 |
Felzartamab | 3.842 | 0.066 | 0.223 | ||
Followup:Felzartamab | 1.915 | 0.162 | 0.347 | ||
Injury-repair induced, day 3 (IRITD3) | Followup | 0.174 | 0.841 | 0.946 | |
Felzartamab | 3.268 | 0.087 | 0.223 | ||
Followup:Felzartamab | 3.125 | 0.056 | 0.223 | ||
Injury-repair induced, day 5 (IRITD5) | Followup | 0.073 | 0.930 | 0.946 | |
Felzartamab | 1.088 | 0.311 | 0.583 | ||
Followup:Felzartamab | 2.586 | 0.089 | 0.223 | ||
Fibrosis classifier (ci>1Prob) | Followup | 0.360 | 0.700 | 0.876 | |
Felzartamab | 3.482 | 0.078 | 0.223 | ||
Followup:Felzartamab | 0.938 | 0.401 | 0.668 | ||
Atrophy classifier (ct>1Prob) | Followup | 0.568 | 0.572 | 0.780 | |
Felzartamab | 3.297 | 0.086 | 0.223 | ||
Followup:Felzartamab | 0.602 | 0.553 | 0.780 | ||
parenchyma-related | Kidney parenchymal (KT1) | Followup | 0.968 | 0.389 | 0.518 |
Felzartamab | 1.628 | 0.218 | 0.518 | ||
Followup:Felzartamab | 0.605 | 0.552 | 0.552 | ||
Kidney parenchymal - no solute carriers (KT2) | Followup | 0.859 | 0.432 | 0.518 | |
Felzartamab | 0.911 | 0.352 | 0.518 | ||
Followup:Felzartamab | 1.038 | 0.365 | 0.518 | ||
Here we
present the effect of treatment (ΔΔ) for each score, along with the
effect of time (Δ) for placebo and felzartmab arms.
flextable_pairwise
Table 1. Effect of felzartamab on molecular scores | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Annotation | Score | Baseline - Week 24 | Week 24 - Week 52 | Baseline - Week 52 | |||||||||
Δ Placebo | Δ Felzartamab | ΔΔ | ΔΔ FDR | Δ Placebo | Δ Felzartamab | ΔΔ | ΔΔ FDR | Δ Placebo | Δ Felzartamab | ΔΔ | ΔΔ FDR | ||
cfDNA | Donor-derived cell-free DNA (dd-cfDNA, cp/mL) | 0.50 (76.12) | -55.00 (30) | -55.50 (53.06) | 2e-04 | -3.50 (69.12) | 14.50 (25.5) | 18.00 (47.31) | 0.027 | -3.00 (51) | -40.50 (45.75) | -37.50 (48.38) | 0.045 |
ABMR-related | ABMR classifier (ABMRProb) | 0.10 (0.37) | -0.57 (0.33) | -0.67 (0.35) | 1e-03 | -0.26 (0.44) | 0.54 (0.45) | 0.79 (0.45) | 1e-03 | -0.15 (0.43) | -0.03 (0.42) | 0.12 (0.43) | 0.709 |
Glomerulitis classifier (g>0Prob) | -0.02 (0.23) | -0.31 (0.35) | -0.29 (0.29) | 3e-03 | -0.19 (0.34) | 0.31 (0.41) | 0.50 (0.37) | 2e-04 | -0.22 (0.3) | 0.00 (0.38) | 0.21 (0.34) | 0.277 | |
Peritubular capillaritis classifier (ptc>0Prob) | -0.12 (0.25) | -0.55 (0.27) | -0.43 (0.26) | 3e-03 | -0.10 (0.37) | 0.37 (0.31) | 0.47 (0.34) | 3e-03 | -0.22 (0.3) | -0.18 (0.35) | 0.04 (0.33) | 0.977 | |
DSA-selective transcripts (DSAST) | -0.09 (0.22) | -0.37 (0.24) | -0.28 (0.23) | 0.026 | -0.16 (0.32) | 0.16 (0.33) | 0.32 (0.33) | 0.026 | -0.24 (0.33) | -0.20 (0.25) | 0.04 (0.29) | 0.988 | |
NK cell burden (NKB) | 0.02 (0.22) | -0.78 (0.49) | -0.80 (0.36) | 0.014 | -0.19 (0.4) | 0.68 (0.59) | 0.87 (0.49) | 0.014 | -0.17 (0.46) | -0.10 (0.62) | 0.07 (0.54) | 0.972 | |
TCMR-related | TCMR classifier (TCMRProb) | -0.01 (0.02) | 0.00 (0) | 0.01 (0.01) | 0.358 | 0.00 (0.02) | 0.00 (0) | 0.00 (0.01) | 0.856 | -0.01 (0.02) | 0.00 (0) | 0.01 (0.01) | 0.358 |
Tubulitis classifier (t>1Prob) | -0.01 (0.06) | -0.02 (0.04) | -0.01 (0.05) | 1.000 | 0.00 (0.04) | 0.01 (0.02) | 0.01 (0.03) | 1.000 | -0.01 (0.05) | -0.01 (0.03) | -0.01 (0.04) | 1.000 | |
Interstitial infiltrate classifier (i>1Prob) | -0.01 (0.04) | -0.02 (0.03) | -0.01 (0.04) | 0.964 | 0.01 (0.04) | 0.01 (0.04) | 0.00 (0.04) | 0.964 | 0.00 (0.05) | -0.01 (0.06) | -0.01 (0.05) | 0.964 | |
Cytotoxic T cell-associated (QCAT) | -0.20 (0.41) | -0.70 (0.46) | -0.50 (0.43) | 0.425 | 0.10 (0.33) | 0.54 (0.57) | 0.44 (0.45) | 0.425 | -0.10 (0.46) | -0.16 (0.41) | -0.06 (0.44) | 0.871 | |
T cell burden (TCB) | -0.44 (0.71) | -0.79 (0.56) | -0.36 (0.64) | 0.825 | 0.28 (0.57) | 0.34 (0.6) | 0.06 (0.58) | 0.905 | -0.16 (0.88) | -0.45 (0.43) | -0.29 (0.65) | 0.825 | |
macrophage-related | Alternatively activated macrophage (AMAT1) | 0.08 (0.2) | -0.27 (0.45) | -0.35 (0.33) | 0.128 | 0.09 (0.23) | 0.11 (0.47) | 0.02 (0.35) | 0.812 | 0.18 (0.19) | -0.16 (0.32) | -0.33 (0.25) | 0.128 |
Constitutive macrophage (QCMAT) | 0.00 (0.15) | -0.24 (0.27) | -0.24 (0.21) | 0.340 | 0.01 (0.12) | 0.09 (0.27) | 0.08 (0.2) | 0.685 | 0.01 (0.11) | -0.15 (0.24) | -0.16 (0.18) | 0.348 | |
injury-related | Injury-repair associated (IRRAT30) | 0.18 (0.52) | 0.05 (0.4) | -0.13 (0.46) | 0.353 | 0.05 (0.66) | -0.11 (0.72) | -0.16 (0.69) | 0.353 | 0.23 (0.6) | -0.06 (0.48) | -0.29 (0.54) | 0.175 |
Injury-repair induced, day 3 (IRITD3) | 0.03 (0.19) | 0.04 (0.12) | 0.01 (0.15) | 0.945 | 0.03 (0.18) | -0.18 (0.18) | -0.20 (0.18) | 0.060 | 0.06 (0.19) | -0.13 (0.14) | -0.19 (0.16) | 0.060 | |
Injury-repair induced, day 5 (IRITD5) | -0.01 (0.17) | 0.04 (0.14) | 0.05 (0.15) | 0.932 | 0.07 (0.2) | -0.17 (0.21) | -0.24 (0.21) | 0.093 | 0.06 (0.18) | -0.13 (0.17) | -0.18 (0.18) | 0.093 | |
Fibrosis classifier (ci>1Prob) | -0.02 (0.29) | 0.13 (0.32) | 0.15 (0.3) | 0.434 | 0.06 (0.25) | -0.08 (0.53) | -0.14 (0.39) | 0.844 | 0.04 (0.25) | 0.05 (0.37) | 0.01 (0.31) | 0.434 | |
Atrophy classifier (ct>1Prob) | 0.01 (0.33) | 0.12 (0.31) | 0.12 (0.32) | 0.587 | 0.09 (0.27) | -0.06 (0.44) | -0.15 (0.36) | 0.587 | 0.10 (0.32) | 0.07 (0.41) | -0.03 (0.36) | 0.587 | |
parenchyma-related | Kidney parenchymal (KT1) | -0.05 (0.22) | -0.15 (0.23) | -0.09 (0.22) | 0.731 | -0.01 (0.11) | 0.10 (0.22) | 0.11 (0.17) | 0.705 | -0.06 (0.25) | -0.04 (0.18) | 0.02 (0.22) | 0.705 |
Kidney parenchymal - no solute carriers (KT2) | -0.09 (0.39) | -0.16 (0.37) | -0.07 (0.38) | 0.577 | -0.09 (0.2) | 0.14 (0.38) | 0.23 (0.29) | 0.484 | -0.18 (0.42) | -0.02 (0.31) | 0.16 (0.37) | 0.577 | |
Grey shading denotes ANOVA interactive effect FDR < 0.05 | |||||||||||||
# EXTRACT LEGEND FOR PLOTS ####
panel_legend <- plots_artANOVA %>%
dplyr::filter(variable == "AMAT1") %>%
pull(plot_violin) %>%
ggpubr::get_legend() %>%
ggpubr::as_ggplot() +
theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))
# MOLECULAR ABMR PANELS ####
panel_violin_abmr <- plots_artANOVA %>%
dplyr::filter(category %in% c("ABMR")) %>%
pull(plot_violin) %>%
wrap_plots(nrow = 1, ncol = 5) &
theme(
axis.text = element_text(size = 10, colour = "black"),
legend.position = "none",
plot.background = element_rect(fill = "grey95", colour = " white")
)
panel_violin_abmr <- panel_violin_abmr %>%
ggarrange(
labels = "B",
font.label = list(size = 25, face = "bold"),
legend = "none"
) %>%
ggpubr::annotate_figure(
top = text_grob("Effect of felzartamab on molecular ABMR activity scores",
face = "bold.italic",
size = 25,
hjust = 1.27
)
)
# MOLECULAR TCMR PANELS ####
panel_violin_tcmr <- plots_artANOVA %>%
dplyr::filter(category %in% c("TCMR")) %>%
pull(plot_violin) %>%
wrap_plots(nrow = 1, ncol = 5) &
theme(
axis.text = element_text(size = 10, colour = "black"),
legend.position = "none",
)
panel_violin_tcmr <- panel_violin_tcmr %>%
ggarrange(
labels = "C",
font.label = list(size = 25, face = "bold"),
legend = "none"
) %>%
ggpubr::annotate_figure(
top = text_grob("Effect of felzartamab on molecular TCMR activity scores",
face = "bold.italic",
size = 25,
hjust = 1.27
)
)
# COMBINED VIOLIN PANELS ####
panels_violin <- ggarrange(
panel_violin_abmr,
panel_violin_tcmr,
nrow = 2,
heights = c(1, 1)
)
ggarrange(
panel_legend,
panels_violin,
nrow = 2,
heights = c(0.125, 1)
)
This analysis
assesses how felzartamab therapy affects gene expression from biopsies
from kidney allografts taken from patients with antibody mediated
rejection. We use the Linear Models for Microarray Data (limma) pipeline
for applying empirical Bayes mediated t-tests on individual genes.
Specific contrasts to assessed the interaction between treatment and
time are specified to estimate the effect of treatemnt on each gene.
References:
• Ritchie, M.E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015).
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# CRAN libraries
library(tidyverse)
library(flextable)
library(officer)
# Bioconductor libraries
library(Biobase)
library(limma)
library(genefilter)
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/data/affymap219.RData")
# Seed for reproducibility
set.seed(42)
# DEFINE THE SET ####
set00 <- set[, set$Patient %nin% c(15, 18)]
• We filter the gene expression data by to remove genes with minimal variance across the population (IQR filtering).
# IQR FILTER THE DATA ####
f1 <- function(x) (IQR(x) > 0.5)
ff <- filterfun(f1)
if (!exists("selected")) {
selected <- genefilter(set00, ff)
}
set01 <- set00[selected, ]
• We then retain a single probeset for each gene based on whichever has the highest mean expression in the population.
# KEEP UNIQUE GENES (keep probe with highest mean expression) ####
mean_exprs_by_probe <- set01 %>%
exprs() %>%
as_tibble(rownames = "AffyID") %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
mutate(mean_exprs = set01 %>%
exprs() %>% rowMeans(), .after = Symb)
genes <- mean_exprs_by_probe %>%
group_by(Symb) %>%
dplyr::slice_max(mean_exprs) %>%
dplyr::filter(Symb != "") %>%
distinct(Symb, .keep_all = TRUE) %>%
pull(AffyID)
set <- set01[featureNames(set01) %in% genes, ]
The
felzartamab trial was a phase II randomize control trial (RCT). While
RCTs are the ‘gold standard’ for establishing causality (i.e., the
effect of treatment), they are not perfect, particularly when the number
of study participants is low. Thus, we pre-screened gene expression data
to remove any genes that differed between placebo and felzartamab arm as
baseline. This was done to reduce the potential for interpreting the
effect of treatment on genes with dissimilar baseline
profiles.
# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set$Felzartamab_Followup %>% droplevels()
# BLOCK DESIGN week24 - baseline ####
design_block01 <- model.matrix(~ 0 + Felzartamab_Followup)
contrast_block_01 <- makeContrasts(
"x = (Felzartamab_FollowupBaseline_Felzartamab) -(Felzartamab_FollowupBaseline_Placebo)",
levels = design_block01
)
# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_block_1 <- limma::lmFit(set, design_block01)
cfit_block_1 <- limma::contrasts.fit(fit_block_1, contrast_block_01)
ebayes_block_1 <- limma::eBayes(cfit_block_1)
tab_block_1 <- limma::topTable(ebayes_block_1, adjust = "fdr", sort.by = "p", number = "all")
# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means_baseline <- fit_block_1 %>%
avearrays() %>%
data.frame() %>%
rownames_to_column("AffyID") %>%
tibble() %>%
mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
dplyr::select(AffyID, contains("Baseline"))
The top
20 genes that differ at baseline (sorted by p-value) between placebo and
felzartamab arms are summarized below.
(Note that all FDR > 0.05. This is a consequence of differential expression among with n = 10 by treatment group across a large selection of genes (5,267 in this case). Nonetheless, genes with p < 0.05 will be excluded from subsequent analyses).
# FORMAT TOPTABLES ####
baseline_deg <- tab_block_1 %>%
as_tibble(rownames = "AffyID") %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID") %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means_baseline, by = "AffyID") %>%
dplyr::select(
AffyID, Symb,
Baseline_Placebo, Baseline_Felzartamab,
t, logFC, P.Value, adj.P.Val,
) %>%
dplyr::rename(p = P.Value, FDR = adj.P.Val)
baseline_deg %>%
dplyr::slice_min(p, n = 20) %>%
mutate_at(vars(t, logFC, FDR), ~ formatC(.x, format = "f", digits = 2)) %>%
mutate_at(vars(p), ~ formatC(.x, format = "f", digits = 5)) %>%
flextable::flextable()
AffyID | Symb | Baseline_Placebo | Baseline_Felzartamab | t | logFC | p | FDR |
|---|---|---|---|---|---|---|---|
11727331_at | SPAG4 | 90 | 50 | -4.23 | -0.85 | 0.00008 | 0.33 |
11763725_a_at | HNRNPU | 366 | 665 | 4.10 | 0.86 | 0.00013 | 0.33 |
11718286_a_at | UBL3 | 74 | 142 | 3.66 | 0.94 | 0.00053 | 0.45 |
11732590_x_at | TRIM5 | 42 | 63 | 3.65 | 0.59 | 0.00055 | 0.45 |
11725671_a_at | IPO7 | 42 | 66 | 3.54 | 0.64 | 0.00079 | 0.45 |
11741431_a_at | KCNJ16 | 203 | 311 | 3.49 | 0.61 | 0.00091 | 0.45 |
11726915_at | SLCO4C1 | 87 | 176 | 3.48 | 1.02 | 0.00096 | 0.45 |
11723109_a_at | GPBP1 | 153 | 263 | 3.46 | 0.78 | 0.00099 | 0.45 |
11739744_x_at | MXI1 | 123 | 191 | 3.45 | 0.64 | 0.00102 | 0.45 |
11720890_a_at | AGL | 39 | 63 | 3.40 | 0.70 | 0.00121 | 0.45 |
11746431_a_at | FBXO11 | 28 | 45 | 3.38 | 0.68 | 0.00127 | 0.45 |
11751841_a_at | TNFRSF17 | 56 | 20 | -3.36 | -1.51 | 0.00134 | 0.45 |
11744367_a_at | ZMYND11 | 30 | 49 | 3.33 | 0.68 | 0.00150 | 0.45 |
11748857_a_at | AK3 | 371 | 578 | 3.32 | 0.64 | 0.00155 | 0.45 |
11727307_x_at | ZNF320 | 46 | 66 | 3.31 | 0.51 | 0.00160 | 0.45 |
11747648_a_at | TCAIM | 96 | 152 | 3.29 | 0.66 | 0.00168 | 0.45 |
11728966_s_at | ZNF28 | 29 | 67 | 3.27 | 1.21 | 0.00181 | 0.45 |
11723387_a_at | PTAR1 | 149 | 261 | 3.25 | 0.81 | 0.00187 | 0.45 |
11739771_a_at | OXR1 | 31 | 51 | 3.24 | 0.73 | 0.00194 | 0.45 |
11730038_at | GNG13 | 44 | 28 | -3.24 | -0.63 | 0.00196 | 0.45 |
We now
carry out differential expression analysis to estimate the treatment
effect on individual genes.
# DEFINE GENES SIMILAR AT BASELINE ####
genes_baseline <- baseline_deg %>%
dplyr::filter(p > 0.05) %>%
pull(AffyID)
# WRANGLE THE MEAN EXPRESSION DATA ####
mean_exprs_by_probe <- set01 %>%
exprs() %>%
as_tibble(rownames = "AffyID") %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
mutate(
mean_exprs = set01 %>%
exprs() %>%
rowMeans(),
.after = Symb
)
genes <- mean_exprs_by_probe %>%
group_by(Symb) %>%
dplyr::slice_max(mean_exprs) %>%
dplyr::filter(Symb != "", AffyID %in% genes_baseline) %>%
distinct(Symb, .keep_all = TRUE) %>%
pull(AffyID)
set02 <- set01[featureNames(set01) %in% genes, ]
We need
to define the model structure for the contrasts we intend to make. Here
we will define 3 major contrast groups; 1) effect of time in placebo
group, 2) effect of time in felzartamab group, and 3) the interactive
effect of time and treatment. We do this for each time window (baseline
to week 24, week 24 to week 52, and baseline to week
52)
# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set02$Felzartamab_Followup %>% droplevels()
cortex <- set02$Cortexprob
# DESIGN ####
design <- model.matrix(~ 0 + Felzartamab_Followup + cortex)
# CONTRAST DESIGN week24 - baseline ####
contrast_interaction_01 <- makeContrasts(
"x = (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 -(Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
contrast_felzartamab_01 <- makeContrasts(
"x = (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
levels = design
)
contrast_placebo_01 <- makeContrasts(
"x = (Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
# CONTRAST DESIGN week52 - week24 ####
contrast_interaction_02 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
levels = design
)
contrast_felzartamab_02 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2",
levels = design
)
contrast_placebo_02 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
levels = design
)
# CONTRAST DESIGN week52 - baseline####
contrast_interaction_03 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
contrast_felzartamab_03 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
levels = design
)
contrast_placebo_03 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
We now
fit the models we’ve specified
# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_interaction_1 <- limma::lmFit(set02, design)
cfit_interaction_1 <- limma::contrasts.fit(fit_interaction_1, contrast_interaction_01)
ebayes_interaction_1 <- limma::eBayes(cfit_interaction_1)
tab_interaction_1 <- limma::topTable(ebayes_interaction_1, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_interaction_1$stdev.unscaled * sqrt(ebayes_interaction_1$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_felzartamab_1 <- limma::lmFit(set02, design)
cfit_felzartamab_1 <- limma::contrasts.fit(fit_felzartamab_1, contrast_felzartamab_01)
ebayes_felzartamab_1 <- limma::eBayes(cfit_felzartamab_1)
tab_felzartamab_1 <- limma::topTable(ebayes_felzartamab_1, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_felzartamab_1$stdev.unscaled * sqrt(ebayes_felzartamab_1$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_placebo_1 <- limma::lmFit(set02, design)
cfit_placebo_1 <- limma::contrasts.fit(fit_placebo_1, contrast_placebo_01)
ebayes_placebo_1 <- limma::eBayes(cfit_placebo_1)
tab_placebo_1 <- limma::topTable(ebayes_placebo_1, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_placebo_1$stdev.unscaled * sqrt(ebayes_placebo_1$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
# FIT BLOCK week52 - week24 LIMMA MODEL ####
fit_interaction_2 <- limma::lmFit(set02, design)
cfit_interaction_2 <- limma::contrasts.fit(fit_interaction_2, contrast_interaction_02)
ebayes_interaction_2 <- limma::eBayes(cfit_interaction_2)
tab_interaction_2 <- limma::topTable(ebayes_interaction_2, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_interaction_2$stdev.unscaled * sqrt(ebayes_interaction_2$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_felzartamab_2 <- limma::lmFit(set02, design)
cfit_felzartamab_2 <- limma::contrasts.fit(fit_felzartamab_2, contrast_felzartamab_02)
ebayes_felzartamab_2 <- limma::eBayes(cfit_felzartamab_2)
tab_felzartamab_2 <- limma::topTable(ebayes_felzartamab_2, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_felzartamab_2$stdev.unscaled * sqrt(ebayes_felzartamab_2$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_placebo_2 <- limma::lmFit(set02, design)
cfit_placebo_2 <- limma::contrasts.fit(fit_placebo_2, contrast_placebo_02)
ebayes_placebo_2 <- limma::eBayes(cfit_placebo_2)
tab_placebo_2 <- limma::topTable(ebayes_placebo_2, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_placebo_2$stdev.unscaled * sqrt(ebayes_placebo_2$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
# FIT BLOCK week52 - baseline LIMMA MODEL ####
fit_interaction_3 <- limma::lmFit(set02, design)
cfit_interaction_3 <- limma::contrasts.fit(fit_interaction_3, contrast_interaction_03)
ebayes_interaction_3 <- limma::eBayes(cfit_interaction_3)
tab_interaction_3 <- limma::topTable(ebayes_interaction_3, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_interaction_3$stdev.unscaled * sqrt(ebayes_interaction_3$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_felzartamab_3 <- limma::lmFit(set02, design)
cfit_felzartamab_3 <- limma::contrasts.fit(fit_felzartamab_3, contrast_felzartamab_03)
ebayes_felzartamab_3 <- limma::eBayes(cfit_felzartamab_3)
tab_felzartamab_3 <- limma::topTable(ebayes_felzartamab_3, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_felzartamab_3$stdev.unscaled * sqrt(ebayes_felzartamab_3$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_placebo_3 <- limma::lmFit(set02, design)
cfit_placebo_3 <- limma::contrasts.fit(fit_placebo_3, contrast_placebo_03)
ebayes_placebo_3 <- limma::eBayes(cfit_placebo_3)
tab_placebo_3 <- limma::topTable(ebayes_placebo_3, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_placebo_3$stdev.unscaled * sqrt(ebayes_placebo_3$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
The
results are now organized for tabulation
# MERGE BLOCKS ####
tab_block_1 <- tab_interaction_1 %>%
left_join(
tab_felzartamab_1 %>%
dplyr::select(AffyID, logFC) %>%
rename(flogFC = logFC),
by = "AffyID"
) %>%
left_join(
tab_placebo_1 %>%
dplyr::select(AffyID, logFC) %>%
rename(plogFC = logFC),
by = "AffyID"
) %>%
relocate(plogFC, flogFC, .before = logFC) %>%
arrange(P.Value)
tab_block_2 <- tab_interaction_2 %>%
left_join(
tab_felzartamab_2 %>%
dplyr::select(AffyID, logFC) %>%
rename(flogFC = logFC),
by = "AffyID"
) %>%
left_join(
tab_placebo_2 %>%
dplyr::select(AffyID, logFC) %>%
rename(plogFC = logFC),
by = "AffyID"
) %>%
relocate(plogFC, flogFC, .before = logFC) %>%
arrange(P.Value)
tab_block_3 <- tab_interaction_3 %>%
left_join(
tab_felzartamab_3 %>%
dplyr::select(AffyID, logFC) %>%
rename(flogFC = logFC),
by = "AffyID"
) %>%
left_join(
tab_placebo_3 %>%
dplyr::select(AffyID, logFC) %>%
rename(plogFC = logFC),
by = "AffyID"
) %>%
relocate(plogFC, flogFC, .before = logFC) %>%
arrange(P.Value)
# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means <- fit_interaction_1 %>%
avearrays() %>%
as_tibble(rownames = "AffyID") %>%
mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
dplyr::select(-contains("Week12"), -any_of(c("cortex")))
# FORMAT TOPTABLES ####
table_interaction_1 <- tab_block_1 %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means, by = "AffyID") %>%
dplyr::select(
AffyID, Symb, Gene, PBT,
t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
all_of(colnames(means))
) %>%
dplyr::rename(
p = P.Value,
FDR = adj.P.Val
)
table_interaction_2 <- tab_block_2 %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means, by = "AffyID") %>%
dplyr::select(
AffyID, Symb, Gene, PBT,
t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
all_of(colnames(means))
) %>%
dplyr::rename(
p = P.Value,
FDR = adj.P.Val
)
table_interaction_3 <- tab_block_3 %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means, by = "AffyID") %>%
dplyr::select(
AffyID, Symb, Gene, PBT,
t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
all_of(colnames(means))
) %>%
dplyr::rename(
p = P.Value,
FDR = adj.P.Val
)
deg <- tibble(
design = c(
"Baseline_vs_Week24",
"Week24_vs_Week52",
"Baseline_vs_Week52"
),
table = list(
table_interaction_1,
table_interaction_2,
table_interaction_3
)
)
# FORMAT TABLES TO MAKE FLEXTABLES ####
table_deg <- deg %>%
expand_grid(direction = c("all", "increased", "decreased")) %>%
relocate(direction, .after = design) %>%
mutate(
gene_tables = pmap(
list(direction, table),
function(direction, table) {
if (direction == "increased") {
table <- table %>%
dplyr::filter(logFC > 0) %>%
arrange(p)
} else if (direction == "decreased") {
table <- table %>%
dplyr::filter(logFC < 0) %>%
arrange(p)
}
table %>%
dplyr::select(-t, -se, -FDR, -contains("AffyID"), -contains("MMDx")) %>%
dplyr::slice(1:20) %>%
mutate(
Gene = Gene %>% str_remove("///.*"),
plogFC = plogFC %>% round(2),
flogFC = flogFC %>% round(2),
logFC = logFC %>% round(2),
p = case_when(
p < 0.0001 ~ p %>% formatC(digits = 0, format = "e"),
TRUE ~ p %>% formatC(digits = 4, format = "f")
)
)
}
)
)
# GLOBAL PARAMETERS FOR FLEXTABLES ####
header1 <- c(
"Gene\nsymbol", "Gene", "PBT",
rep("Differential expression", 4),
rep("Mean expression by group", 6)
)
header2 <- c(
"Gene\nsymbol", "Gene", "PBT",
"\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
"\u394\u394\nlogFC", "\u394\u394\nP",
rep("Placebo", 3), rep("Felzartamab", 3)
)
header3 <- c(
"Gene\nsymbol", "Gene", "PBT",
"\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
"\u394\u394\nlogFC", "\u394\u394\nP",
"Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)", "Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)"
)
# MAKE FORMATTED FLEXTABLES ####
flextables <- table_deg %>%
mutate(
flextable = pmap(
list(design, gene_tables),
function(design, gene_tables) {
if (design == "Baseline_vs_Week24") {
title <- paste("Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
} else if (design == "Week24_vs_Week52") {
title <- paste("Table i. Top 20 differentially expressed genes between week24 and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
} else if (design == "Baseline_vs_Week52") {
title <- paste("Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
}
gene_tables %>%
arrange(logFC) %>%
flextable::flextable() %>%
flextable::delete_part("header") %>%
flextable::add_header_row(top = TRUE, values = header3) %>%
flextable::add_header_row(top = TRUE, values = header2) %>%
flextable::add_header_row(top = TRUE, values = header1) %>%
flextable::add_header_row(top = TRUE, values = rep(title, ncol_keys(.))) %>%
flextable::merge_v(part = "header") %>%
flextable::merge_h(part = "header") %>%
flextable::border_remove() %>%
flextable::border(part = "header", border = officer::fp_border()) %>%
flextable::border(part = "body", border = officer::fp_border()) %>%
flextable::align(align = "center") %>%
flextable::align(align = "center", part = "header") %>%
flextable::font(fontname = "Arial", part = "all") %>%
flextable::fontsize(size = 7, part = "all") %>%
flextable::fontsize(size = 7, part = "footer") %>%
flextable::fontsize(i = 1, size = 12, part = "header") %>%
flextable::bold(part = "header") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::padding(padding = 0, part = "all") %>%
flextable::autofit()
}
)
)
Genes
suppressed by felzartamab in week 24 compared to baseline
(Note
that all ΔΔFDR > 0.05. This is a consequence of differential
expression among with n = 10 by treatment group across a large selection
of genes (5,267 in this case). Thus, treatment effects on individual
genes should not be interpreted. Instead, geneset enrichment analysis
(GSEA) will be used to assess bulk trends).
flextables %>%
dplyr::filter(design == "Baseline_vs_Week24", direction == "decreased") %>%
pull(flextable) %>%
pluck(1)
Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value) | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
Gene | Gene | PBT | Differential expression | Mean expression by group | ||||||||
Δ | Δ | ΔΔ | ΔΔ | Placebo | Felzartamab | |||||||
Baseline | Week24 | Week52 | Baseline | Week24 | Week52 | |||||||
CXCL11 | chemokine (C-X-C motif) ligand 11 | ABMR-RAT,GRIT3,RAT,Rej-RAT | -0.16 | -1.01 | -0.85 | 0.0176 | 494 | 395 | 372 | 617 | 153 | 317 |
CXCL9 | chemokine (C-X-C motif) ligand 9 | ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT | -0.14 | -0.94 | -0.79 | 0.0290 | 1,344 | 1,102 | 1,384 | 1,684 | 459 | 1,036 |
CXCL10 | chemokine (C-X-C motif) ligand 10 | ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT | -0.21 | -0.98 | -0.77 | 0.0228 | 844 | 632 | 720 | 1,063 | 274 | 615 |
IDO1 | indoleamine 2,3-dioxygenase 1 | ABMR-RAT,GRIT3,RAT,Rej-RAT | -0.09 | -0.77 | -0.68 | 0.0103 | 338 | 298 | 296 | 409 | 140 | 251 |
LYZ | lysozyme | IRITD5,QCMAT | -0.01 | -0.67 | -0.66 | 0.0386 | 1,676 | 1,658 | 1,884 | 2,244 | 884 | 1,184 |
FCGR3A | Fc fragment of IgG, low affinity IIIa, receptor (CD16a) | ABMR-RAT,GRIT2,IRRAT950,RAT,Rej-RAT | 0.06 | -0.58 | -0.64 | 0.0167 | 550 | 602 | 613 | 571 | 256 | 352 |
GZMB | granzyme B | ABMR-RAT,QCAT,RAT,Rej-RAT | -0.01 | -0.62 | -0.61 | 0.0019 | 149 | 147 | 130 | 143 | 61 | 114 |
GBP1 | guanylate binding protein 1, interferon-inducible | ABMR-RAT,GRIT3,RAT,Rej-RAT | -0.03 | -0.61 | -0.59 | 0.0096 | 678 | 655 | 643 | 909 | 389 | 556 |
GNLY | granulysin | ABMR-RAT,DSAST,QCAT,RAT,Rej-RAT | -0.18 | -0.75 | -0.57 | 0.0285 | 192 | 149 | 134 | 165 | 58 | 100 |
SH2D1B | SH2 domain containing 1B | ABMR-RAT,DSAST,NKB,RAT | 0.07 | -0.44 | -0.51 | 0.0010 | 38 | 42 | 34 | 40 | 22 | 33 |
FCGR1A | Fc fragment of IgG, high affinity Ia, receptor (CD64) | GRIT2,IRRAT950,RAT,TCMR-RAT | 0.09 | -0.41 | -0.50 | 0.0304 | 131 | 148 | 161 | 152 | 86 | 86 |
KLRD1 | killer cell lectin-like receptor subfamily D, member 1 | ABMR-RAT,RAT,Rej-RAT | -0.06 | -0.49 | -0.43 | 0.0295 | 110 | 101 | 96 | 126 | 63 | 102 |
EPB41L3 | erythrocyte membrane protein band 4.1-like 3 | 0.22 | -0.20 | -0.42 | 0.0388 | 209 | 286 | 258 | 339 | 258 | 283 | |
LOC339059 | uncharacterized LOC339059 | 0.26 | -0.15 | -0.41 | 0.0203 | 51 | 73 | 55 | 66 | 53 | 59 | |
APOL1 | apolipoprotein L1 | ABMR-RAT,GRIT3,RAT,Rej-RAT | 0.07 | -0.32 | -0.39 | 0.0205 | 626 | 687 | 744 | 754 | 483 | 636 |
LYPD5 | LY6/PLAUR domain containing 5 | ABMR-RAT,RAT,Rej-RAT | 0.04 | -0.33 | -0.37 | 0.0255 | 27 | 28 | 27 | 30 | 19 | 26 |
TRAV12-2 | T cell receptor alpha variable 12-2 | 0.25 | -0.10 | -0.35 | 0.0173 | 39 | 55 | 39 | 49 | 43 | 43 | |
IL18RAP | interleukin 18 receptor accessory protein | 0.02 | -0.32 | -0.34 | 0.0233 | 38 | 39 | 41 | 42 | 27 | 33 | |
PRF1 | perforin 1 (pore forming protein) | ABMR-RAT,QCAT,RAT,Rej-RAT | -0.02 | -0.35 | -0.33 | 0.0382 | 152 | 148 | 137 | 145 | 89 | 123 |
PMM2 | phosphomannomutase 2 | 0.20 | -0.07 | -0.27 | 0.0380 | 38 | 50 | 41 | 46 | 42 | 49 | |
Genes
suppressed by felzartamab in week 52 compared to baseline
(Note
that all ΔΔFDR > 0.05. This is a consequence of differential
expression among with n = 10 by treatment group across a large selection
of genes (5,267 in this case). Thus, treatment effects on individual
genes should not be interpreted. Instead, geneset enrichment analysis
(GSEA) will be used to assess bulk trends).
flextables %>%
dplyr::filter(design == "Baseline_vs_Week52", direction == "decreased") %>%
pull(flextable) %>%
pluck(1)
Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value) | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
Gene | Gene | PBT | Differential expression | Mean expression by group | ||||||||
Δ | Δ | ΔΔ | ΔΔ | Placebo | Felzartamab | |||||||
Baseline | Week24 | Week52 | Baseline | Week24 | Week52 | |||||||
HNRNPK | heterogeneous nuclear ribonucleoprotein K | 0.28 | -0.38 | -0.66 | 0.0055 | 818 | 1,044 | 1,203 | 1,096 | 1,059 | 649 | |
COL1A1 | collagen, type I, alpha 1 | COLA1356,FICOL,IRITD5,LivGST_UP | 0.38 | -0.26 | -0.64 | 0.0075 | 310 | 333 | 526 | 329 | 320 | 230 |
NFIX | nuclear factor I/X (CCAAT-binding transcription factor) | CT1 | 0.42 | -0.15 | -0.57 | 0.0012 | 162 | 179 | 289 | 221 | 167 | 180 |
RHOA | ras homolog family member A | IRRAT950 | 0.28 | -0.29 | -0.57 | 0.0049 | 691 | 803 | 1,017 | 923 | 809 | 620 |
COL1A2 | collagen, type I, alpha 2 | FICOL,IRITD5,LivGST_UP | 0.26 | -0.31 | -0.57 | 0.0118 | 1,238 | 1,152 | 1,769 | 1,534 | 1,194 | 997 |
ITGB3 | integrin beta 3 | IRRAT30,IRRAT950 | 0.34 | -0.21 | -0.55 | 0.0100 | 105 | 137 | 169 | 110 | 103 | 82 |
LOC100129518 | uncharacterized LOC100129518 | GRIT3,IRRAT950 | 0.25 | -0.22 | -0.47 | 0.0050 | 429 | 546 | 605 | 527 | 461 | 390 |
ACTB | actin, beta | IRRAT950 | 0.15 | -0.31 | -0.47 | 0.0077 | 6,451 | 7,352 | 7,993 | 7,274 | 5,662 | 4,720 |
SPARC | secreted protein, acidic, cysteine-rich (osteonectin) | IRITD3 | 0.16 | -0.30 | -0.46 | 0.0100 | 2,411 | 2,532 | 3,011 | 2,343 | 2,094 | 1,550 |
DNAJA4 | DnaJ (Hsp40) homolog, subfamily A, member 4 | IRRAT950 | 0.27 | -0.14 | -0.42 | 0.0095 | 49 | 68 | 72 | 68 | 61 | 56 |
SLFN5 | schlafen family member 5 | 0.19 | -0.20 | -0.40 | 0.0050 | 39 | 44 | 51 | 49 | 38 | 37 | |
CYR61 | cysteine-rich, angiogenic inducer, 61 | IRITD3,IRRAT950 | 0.14 | -0.26 | -0.40 | 0.0074 | 271 | 268 | 330 | 338 | 234 | 236 |
SDCBP | syndecan binding protein | 0.24 | -0.14 | -0.38 | 0.0120 | 1,135 | 1,404 | 1,582 | 1,478 | 1,358 | 1,224 | |
SBSPON | somatomedin B and thrombospondin type 1 domain containing | 0.31 | -0.08 | -0.38 | 0.0123 | 217 | 250 | 332 | 188 | 249 | 169 | |
PAFAH1B2 | platelet-activating factor acetylhydrolase 1b, catalytic subunit 2 (30kDa) | 0.23 | -0.13 | -0.36 | 0.0099 | 309 | 372 | 423 | 386 | 396 | 323 | |
TMEM2 | transmembrane protein 2 | 0.10 | -0.26 | -0.36 | 0.0118 | 270 | 307 | 311 | 318 | 259 | 222 | |
CALM2 | calmodulin 2 (phosphorylase kinase, delta) | CT1,IRITD3 | 0.12 | -0.25 | -0.36 | 0.0138 | 1,304 | 1,427 | 1,538 | 1,623 | 1,342 | 1,155 |
WNK1 | WNK lysine deficient protein kinase 1 | 0.26 | -0.08 | -0.34 | 0.0104 | 703 | 884 | 1,012 | 894 | 996 | 799 | |
CEBPB | CCAAT/enhancer binding protein (C/EBP), beta | IRITD3 | 0.11 | -0.22 | -0.33 | 0.0123 | 419 | 452 | 488 | 417 | 333 | 307 |
STAT6 | signal transducer and activator of transcription 6, interleukin-4 induced | 0.09 | -0.20 | -0.30 | 0.0144 | 671 | 710 | 765 | 746 | 616 | 563 | |